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1.  Introduction 


A  central  challenge  engine  manufactures  face  is  meeting  the  power  requirements 
with  ever  higher  combustion  efficiency.  The  US  Army’s  need  for  faster,  leaner, 
more  powerful  engines  is  motivating  the  power  and  energy  industry  to  produce 
technologies  that  can  meet  these  requirements.  Novel  engine  designs  must  provide 
fuel  flexibility,  improved  combustion  efficiencies,  and  power  densities  to  ensure 
battlefield  dominance  for  air  and  terrestrial  operations.  For  Army  operations, 
reducing  the  battlefield  fuel  consumption  (and  payload)  can  have  a  significant 
impact  on  improving  the  military  operational  capability,  enhancing  force 
projection,  and  contributing  to  mission  success. 

In  addressing  these  challenges  it  is  important  to  understand  the  quality  of  the 
fuel-air-mixture  formation  process  as  it  directly  affects  the  combustion 
performance  (Bravo  et  al.  2014b,  2015a,  2015b,  2016a;  Ma  et  al.  2014).  Mixture 
formation  starts  with  the  injection  of  the  compressed  liquid  fuel,  which  mixes  and 
vaporizes  with  the  surrounding  oxidizer  leading  to  ignition  and  sustained 
combustion  (Ray  et  al.  2016).  It  is  important  to  understand  the  behavior  for  the 
liquid  length  in  engine  sprays  in  order  to  develop  mixture-formation  strategies  and 
avoid  adverse  outcomes  such  as  spray  hitting  the  piston  wall  and  causing  increased 
emissions.  To  deal  with  this  trade-off,  a  number  of  spray  models  and  approaches 
have  been  developed.  This  report  compares  2  such  models:  a  reduced-order  model 
(ROM)  developed  from  the  work  of  Naber  and  Siebers  (Naber  and  Siebers  1996; 
Siebers  1999)  and  a  full-order  model  (FOM)  based  on  3-D-computational  fluid 
dynamics  (CFD)  simulations.  The  FOM  is  used  as  an  initial  benchmark,  as  it  has 
several  advantages  including  its  formulation  based  on  fluid  governing  equations, 
multidimensionality,  and  numerical  accuracy.  As  usual  with  such  models,  the  main 
disadvantage  of  the  FOM  simulation  is  its  steep  computational  cost,  long  simulation 
time  (order  of  days),  and  execution  on  high-performance-computing  (HPC) 
platforms  that  may  not  be  readily  available  in  the  field.  By  contrast,  a  ROM  takes 
fewer  inputs  and  can  run  on  a  standard  desktop  or  laptop  equipped  with  a  more 
common  programming  language  (FORTRAN,  C++,  etc.).  The  main  disadvantage 
of  the  ROM  is  that  it  is  limited  to  lower  dimensions  and  steady-state  conditions  and 
introduces  assumptions  (e.g.,  mixing-limited  vaporization  that  limits  its 
applicability). 

The  original  Siebers  model  (Siebers  1999)  introduced  a  mixing-limited  scaling  law 
(zero-dimensional)  for  the  steady-state  calculation  of  an  evaporating  liquid  spray 
plume.  Siebers’  model  combined  measured  input  functions  (mainly  from 
curve-fitting  American  Petroleum  Institute  [API]  enthalpy  and  compressibility 
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data)  with  empirical  parameters  to  observed  liquid  lengths.  This  method  found  a 
good  fit  for  most  points  except  in  the  low-density  range  where  the  law 
overestimated  the  liquid  length.  In  earlier  work,  Siebers  and  Naber  had  introduced 
a  penetration  correlation  that  estimates  the  transient  spray  state  as  non  vaporizing. 
A  paper  by  Schihl  et  al.  (Schihl  et  al.  2006)  later  revisited  and  extended  the  model 
finding  good  agreement  with  experiment  in  the  ranges  of  interest.*  To  expand  the 
model’s  scope,  Schihl  et  al.  also  introduced  several  schemes  for  analyzing  more 
realistic  fuel  blends,  all  based  on  representing  the  fuel  as  a  surrogate  mixture  of 
pure  components.  A  series  of  papers  by  Kurvers  and  Luijten  (Kurvers  2009; 
Kurvers  and  Luijten  2010;  Luijten  and  Kurvers  2010)  examined  the  effect  of  real- 
gas  corrections  to  the  scaling  law.  They  found  some  improvement  when  accounting 
for  real-gas  effects;  however,  the  results  continued  to  marginally  overestimate  the 
penetration  lengths  at  low  densities.  In  summary,  these  studies  suggest  the 
assumptions  behind  the  model  are  largely  valid  at  most  conditions  relevant  to  diesel 
engines,  although  the  models  struggle  at  densities  below  around  14  kg/m3. 

This  report  describes  an  implementation  of  the  Siebers  steady-state  scaling  law 
along  with  the  Naber-Siebers  penetration  correlation  for  transient  effects.  The 
cases  studied  include  an  evaporating  single-plume  spray  and  a  single-cylinder 
moving  piston  case  near  top  dead  center  (TDC)  at  diesel-engine  conditions.  The 
results  are  compared  to  laboratory  measurements  and  to  FOM  simulations.  The 
ROM  is  coded  in  MATLAB,  and  the  FOM  is  modeled  using  the  CONVERGE 
commercial  CFD  package.  The  FOM  requires  the  use  of  a  HPC  platform;  in  this 
case,  the  Garnet  system  at  the  Army  Corps  of  Engineers’  Engineer  Research  and 
Development  Center.  The  ROM  provides  a  real-time  engineering  analytical  tool  for 
liquid-length  scaling  that  may  be  used  toward  optimizing  engine  performance. 

2.  ROM  Methodology 

The  ROM  consists  of  2  segments.  The  first  is  a  steady-state  solution  at 
thermodynamic  equilibrium  based  on  Siebers’  method  (Siebers  1999).  Siebers 
derived  his  model,  which  he  termed  a  scaling  law,  by  combining  conservation  of 
mass,  momentum,  and  energy  in  a  fuel  spray  with  phase-transition  thermodynamics 
at  the  liquid  length.  Implementing  this  law  requires  a  number  of  other  inputs.  These 
may  be  either  empirical  or  based  on  theoretical  concepts.  The  model  presented  here 
largely  follows  the  empirical  approach  of  Siebers  (1999)  and  explored  in  detail  by 
Schihl  et  al.  (2006).  Figure  1  depicts  the  concept  of  Siebers’  framework. 


*  It  should  be  noted  Schihl  et  al.  (2006)  contains  what  seems  to  be  a  typographical  error  in  its  replication 
of  Siebers’  scaling  law.  The  numerator  of  its  Equation  2  contains  a  factor  of  Pa  where  Ps  should  be  used 
instead. 
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A  second  segment  determines  transient  behavior  between  injection  and 
equilibrium.  It  is  based  on  a  penetration  correlation  derived  by  Naber  and  Siebers 
(1996),  drawn  largely  on  the  same  conceptual  framework  as  the  steady-state 
solution  but  without  the  imposition  of  thermodynamic  conditions.  The  flow  chart 
in  Fig.  2  shows  a  schematic  of  how  the  ROM  is  implemented. 

2.1  Sieber's  Theoretical  Framework 

The  Siebers  ROM  is  conceptually  based  on  a  mixing-limited  conical  spray,  which 
can  be  visualized  as  a  continuous  series  of  zero-depth  disks  of  fuel-air  blend  that 
propagate  forward  from  the  nozzle,  expanding  and  entraining  air  as  they  go.  Each 
disk  is  gradually  vaporized  until  the  time  it  reaches  a  distance  L,  where  all  fuel  is 
vaporized.  This  distance  is  referred  to  as  the  steady-state  liquid  length. 


Entrained  gas  (Pa ,  Ta ,  pa) 


complete  (x=L) 

Fig.  1  Schematic  of  theoretical  conceptual  framework  by  Siebers  (1999)  in  research 
sponsored  by  the  US  Department  of  Energy 

The  formulation  begins  with  mass  conservation 

rhf  =  rhf0  =  pfA0U0  ,  (1) 

ma  =  paA  U  (2) 

and  momentum  conservation 

rhf0U0  =  rhfU  +  maU  .  (3) 

Here  A  is  the  area  of  a  disk,  U  its  velocity,  and  m  the  mass  flow  rate  of  the  specified 
component  (/  for  fuel,  a  for  entrained  air).  Densities  are  constant.  Other  quantities 
marked  with  a  zero  are  values  at  the  nozzle.  All  others  are  implicitly  dependent  on 
the  distance  x  from  the  nozzle. 
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Input  file  with  chamber 
conditions  and 
component  mixtures 

A 

Load  component  data  from  subfiles 
and  compute  operating  conditions 


Ca  Ic Li  late  s  p  ray  a  ng  I  e , 
steady-state  liquid  length 


A 

L 


Return  L(t)  forgiven  t 


Fig.  2  Simplified  schematic  flow  of  ROM  calculation 
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Using  the  conical  geometry  and  empirical  values  provided  by  Siebers  (1999),  one 
can  eliminate  the  velocities  to  solve  for  x 


x  =  0.62 


deff 

tana/2 


i  , 


(4) 


with 


g  _  rhfjx) 
ma(x) 


(5) 


The  liquid  length  L  is  just  the  x  value  associated  with  B(L),  but  since  the  formula 
is  self-referencing  and  too  difficult  to  solve,  this  value  must  be  determined  in  turn 
by  iteration.  Siebers  (1999)  uses  thermodynamic  properties  of  the  evaporating 
spray  to  recast  B  in  2  further  ways. 

First  from  the  enthalpy  conservation  between  the  initial  final  states  of  an  air-fuel 
disk, 


rhfhf(jf,pa)  +  rhaha(Ta,Pa )  =  mfhf(Ts )  +  rhaha(Ts,Pa  -  Ps )  6) 

with  specific  enthalpies  ha  and  hf  for  the  air  and  gas  respectively.  Rearranging  and 
assuming  the  fuel-air  mix  reaches  thermal  equilibrium  ( Ts ,  Ps)  at  the  liquid  length 
gives 


_  raf  hq(Tg,Pa)  hgtTsJji  Ps) 

liT-a  hf(Ts)  ~  hf{Tf,Pa) 

In  this  context  B  is  known  as  the  evaporation  constant. 
Second  from  ideal  gas  law 


P  =  Z -m- MW  ■  RT 

with  compressibility  Z  and  molecular  weight  MW.  Since  the  fuel  and  air  are  both 
vapor  at  the  liquid  length,  rearranging  gives  another  form  of  B: 

B=^/=  Za(Ts,Pg-Ps)  Ps  MWf 

ma  Zf(Ts,Ps )  Pa-Ps  MWa  '  1 

Determining  the  evaporation  constant  is  now  equivalent  to  finding  a  Ts  that  satisfies 

hq (Tq)  -  ftq(Ts)  _  Za(Ts)  P5  MW f 

hf(Ts)-hf(Tf)  Zf(Ts )  Pa-Ps  MWa  ‘  ^ 
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The  model  is  completed  by  input  functions  for  the  enthalpies,  compressibility,  and 
saturation  pressure.  The  next  section  will  detail  how  this  evaporation  condition  is 
discovered. 


2.2  Solution  Method 

The  2  representations  of  B  can  be  reformulated  as  2  functions  of  Ts: 

n  /-y  \  _  ^g(^g)  ~  hg(Ts) 

s)  ~  hf(Ts)-hf{Tf)  ’ 

B  (T  )  =  ?a P’s)  .  A(A)  . 

S)  Zf(Ts )  Pa-Ps(W  MWa 

By  comparing  them,  the  evaporation  constant  can  be  determined  for  a  fuel  species 
via  a  root-finding  algorithm  that  adjusts  Ts  until  these  functions  are  within  a 
specified  tolerance  of  each  other: 


g  & !h_ 
B 


<  £ 


(10) 


with  the  average  B  =  ^  (Bg  +  Blr)  taken  as  the  final  value. 

The  MATLAB  algorithm  in  this  study  does  this  by  guessing  a  series  of  values  for 
Ts  and  updating  these  guesses  up  or  down  depending  on  the  difference  found  in 
Eq.  10.  This  algorithm  is  similar  to  a  binary  search  on  an  array,  except  here  the 
search  space  is  the  real  line.  As  in  a  binary  search,  convergence  occurs  in 
logarithmic  time,  as  shown  in  Appendix  A.  Appendix  B  has  the  MATLAB  code  for 
ROM. 


The  input  values  and  tables  of  species-specific  parameters  are  fed  into  files  that  run 
function  files.  The  functions  used  are  empirical  formulas  (which  follow,  with 
references).  Most  of  these  functions  require  further  species-specific  parameters, 
which  are  given  for  the  3  pure  fuels  and  4  ambient  gases  the  ROM  currently  uses 
(see  Tables  1  and  2).  Values  of  h  are  in  joules  per  gram;  all  others  are  in  SI  units 
except  where  specified. 

Ambient  enthalpies — NASA  formula  (McBride  et  al.  1993): 

ha(T  -  298.15  K)  =—  (aT  +  B  —  +  C  —  +  D  —  +  E  —  +  F  -  g)  .  (11) 

av  J  MWa  \  2  3  4  5  /  v  ' 
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Table  1  Species-specific  parameters  for  4  ambient  gases  the  ROM  uses 


Gas 

Range 

A 

BTO3 

CIO7 

DT01# 

ETO14 

FTO'3 

GIO'4 

Nitrogen 

T  <  1000  K 

3.531 

-0.1237 

-5.030 

24.35 

-140.9 

-1.047 

0.000 

T  >  1000  K 

2.953 

1.397 

-4.926 

0.7860 

-0.4608 

-0.9239 

0.000 

Oxygen 

T  <  1000  K 

3.782 

-2.997 

98.47 

-96.81 

324.4 

-1.064 

0.000 

T  >  1000  K 

3.661 

0.6564 

-1.411 

0.2058 

-0.1299 

-1.216 

0.000 

Carbon 

T  <  1000  K 

2.357 

8.985 

-71.23 

24.59 

-14.37 

-48.37 

-4.733 

dioxide 

T  >  1000  K 

4.637 

2.741 

-9.958 

1.604 

-0.9161 

-49.02 

-4.733 

Water 

T  <  1000  K 

4.199 

-2.037 

65.20 

-54.88 

177.2 

-30.29 

-2.908 

vapor 

T  >  1000  K 

2.677 

2.973 

-7.738 

0.9443 

-0.4269 

-29.89 

-2.908 

Saturation  pressure — API  handbook  (API  1997),  Formula  5A1.1*: 

\nPs=  A-  BTf1  -C\nTs  +  DTS2  +  ET~ 2  .  (12) 

Table  2  Species-specific  parameters  for  3  pure  fuels  the  ROM  uses 


Fuel 

A 

B 

c 

DIO6 

E 

Dodecane 

170.27 

25990 

20.822 

2.9759 

645190 

Tetradecane 

130.78 

20072 

15.743 

2.531 

-1008900 

Cetane 

174.2 

28534 

21.09 

2.5228 

88111 

Liquid  density — API  handbook  (API  1997),  Formula  6A2.13-1: 

J_  —  7  1+(1— T //7c)2/7 

r\  ~  d  L ra 


Liquid-fuel  enthalpy — Schihl  et  al.  (2006)  curve  fit  of  API  data: 


h,(T,P)  =  { 


Prs  <  0.2: 
Prs  >  0.2: 


ATrs  -  B 

— CTrs2  +  DTrs  —  E 


Vaporized-fuel  enthalpy — Schihl  et  al.  (2006)  curve  fit  of  API  data: 

hf(Pa)=APf  +  B  . 

rc 

Table  3  lists  Eq.  14’s  model  coefficients  for  the  3  pure  fuels: 

Table  3  Model  coefficients  for  Eq.  14 


(13) 


(14) 


(15) 


Fuel 

A 

B 

Dodecane 

1.277 

316.3 

Tetradecane 

1.1718 

319.2 

Cetane 

1.0429 

320.95 

*  Method  API  5A1.1  takes  temperature  input  in  Rankine  and  outputs  pressure  in  psi  absolute. 
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Fuel-vapor  compressibility — Schihl  et  al.  (2006)  curve  fit  of  API  data: 


Zf  =  —A 


Table  4  lists  Eq.  16’s  model  coefficients  for  the  3  pure  fuels: 


Table  4  Model  coefficients  for  Eq.  16 


Fuel 

A 

B 

c 

D 

Dodecane 

16.85 

36.104 

26.425 

7.5406 

Tetradecane 

17.924 

36.143 

24.71 

6.6857 

Cetane 

16.587 

34.594 

24.531 

6.869 

Spray  angle — Siebers  (1999)  curve  fit  of  measurements*: 


tan- =  0.264 

2 


0.0043 


(16) 


(17) 


Finally,  for  fuel  blends  the  liquid  length  is  determined  through  a  mean  evaporation 
constant  (MEC)  scheme  '  suggested  by  Schihl  et  al.  (2006).  The  properties  are 
calculated  for  each  pure  component  of  the  blend  according  to  Eqs.  1-7.  This  gives 
an  evaporation  constant  B*  for  each  component.  A  weighted  average  gives 

B=  Yi=ixi  Bi  (18) 

where  xt  is  the  mass  fraction  of  each  component.  This  result  is  treated  as  the 
evaporation  constant  for  the  entire  mixture  and  input  into  the  scaling  law  to  return 
the  blend’s  liquid  length. 


2.3  Transient  Modeling 

The  scaling-law  methodology  is  capable  only  of  estimating  the  liquid  length  at 
steady  state.  The  Siebers  framework  only  attempts  to  capture  the  behavior  of  the 
spray  at  steady  state.  To  account  for  transient  behavior,  the  ROM  is  modified 
according  to  the  work  presented  by  Naber  and  Siebers  (1996).  In  their  paper,  Naber 
and  Siebers  use  the  same  conceptual  model  as  in  this  report’s  Section  2.1  to 
introduce  a  correlation  for  spray  growth  over  time.  The  major  new  assumption  of 
this  transient  model  is  to  ignore  vaporization  and  substitute  spray  length  for  liquid 
penetration.  As  a  result,  the  transient  model  relies  entirely  on  the  kinematics  of  the 


*  Strictly  speaking,  the  factor  of  0.264  varies  by  nozzle,  but  this  parameter  showed  no  clear  empirical 
relationship  to  diameter.  For  the  nozzles  studied  by  Siebers  (1999)  this  value  did  not  vary  much,  so  we  use 
the  mean  of  those  nozzles  as  a  default. 

1  Among  other  schemes  tested  by  Schihl  et  al.  (2006),  a  simple  weighted  average  of  the  component  liquid 
lengths  was  also  found  to  be  adequate,  but  Schihl  et  al.  found  the  MEC  scheme  to  be  slightly  more  accurate. 
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spray.  The  resulting  transient  length  is  therefore  only  approximate,  but  this  study 
shows  it  can  be  suitably  modified.  The  derivation  is  otherwise  similar  to  that  of 
steady  state  in  Siebers  (1999).  Naber  and  Siebers  (1996)  first  define  a  length  scale: 


L+  = 


deff 

0.66  tana/2 


(19) 


This  length  scale  also  appears  implicitly  in  the  steady-state  scaling  law.  Time  and 
the  time-dependent  liquid  length  can  be  related  by  reformulating  them  as 

L(t)  =  L+  ■  Lit)  (20) 


and 


t=7T-t  (21) 

uf 

with  normalized  variables  L  and  t.  Here,  Uf  is  the  injection  velocity  of  the  fuel  as 
determined  by  the  Bernoulli  law 


Vf  =  <22> 

given  fuel-injection  pressure  Pf  and  an  empirical- velocity  coefficient  Cv.  The 
Bernoulli  law  applies  to  the  fluid  velocity,  not  the  liquid-length  growth,  so  Cv 
should  be  close  to  unity.  Upon  analysis  of  the  validation  data,  we  found  that 
observed  liquid  lengths  grew  3-5  times  slower  than  the  injection  velocity. 
Adjusting  Cv  amounts  to  modifying  the  velocity  of  the  spray,  which  can  be  done 
until  it  matches  the  empirical  value  of  the  liquid.  Nozzles  examined  in  Siebers’ 
paper  (1999)  had  velocity  coefficients  between  0.8  and  1  for  the  entire  spray.  The 
value  we  use  for  the  liquid  portion  is  0.31,  which  most  closely  approximates  the 
observed  behavior  of  the  liquid  length. 

Applying  the  kinematics  of  the  spray  using  Eqs.  1-3,  one  may  achieve  a  correlation 
between  Eqs.  20  and  21: 

t  =  h-  +  ^Vl  +  16L2  +  ^ln  (4 L  +  V 1  +  16L2)  .  (23) 

This  correlation  is  too  difficult  to  invert,  so  it  must  be  solved  numerically.  Our 
implementation  closely  resembles  the  iterative  algorithm  described  in  the  previous 
section  to  determine  the  evaporation  constant.  Given  a  time  t,  a  guess  is  made  and 
updated  for  L  until  Eq.  23  evaluates  to  a  time  sufficiently  close  to  t. 
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As  noted  by  Naber  and  Siebers  (1996),  the  transient  correlation  is  unbounded,  so  it 
will  not  accurately  demonstrate  steady-state  behavior.  This  is  corrected  with  a 
piecewise  function  in  time.  First,  the  steady-state  liquid  length  is  determined 
according  to  the  Siebers  scaling  law.  The  time  required  to  reach  this  distance  is 
calculated  as  tL,  according  to  Eq.  23.  For  times  before  tL,  liquid  lengths  are 
determined  according  to  the  penetration  correlation.  For  times  after,  the  result  from 
the  scaling  law  is  used. 

3.  FOM  Methodology 

The  CONVERGE  3D-CFD  solver,  developed  by  Convergent  Science,  Inc.,  has 
been  adopted  in  this  study  to  perform  detailed  spray  simulations  at  realistic 
engine-operating  conditions.  CONVERGE  is  a  compressible  Navier-Stokes  solver. 
It  is  based  on  a  first-order  predictor-time-integration  scheme  and  a  choice  of  second 
or  higher  order  finite-volume  schemes  for  spatial  discretization.  It  features  a 
nonstaggered,  collocated  computation  grid  framework  using  a  Rhie-Chow 
interpolation  technique  to  avoid  spurious  oscillations.  An  efficient  geometric 
multigrid  treatment  is  used  to  solve  the  pressure  equation,  and  a  parallel  computing 
implementation  is  based  on  implementations  of  either  OpenMP  or  Message  Passing 
Interface  protocols.  It  provides  the  option  of  increasing  resolution  locally  through 
static  fixed-grid  embedding  and  dynamically  through  Adaptive  Mesh  Refinement 
activated  through  user-specified  criteria.  Additionally,  it  uses  state-of-the-art 
Eulerian-Lagrangian  spray  models  to  describe  particle  transport  such  as  liquid-fuel 
injection  (Bravo  et  al.  2016,  2017). 

In  this  study,  the  gas  phase  is  described  using  the  Favre-Averaged  Navier-Stokes 
equations,  and  the  renormalization  group’s  (RNG’s)  k  —  e  model  for  low- 
resolution  turbulence  studies  has  been  adopted.  The  compressible  system  of 
transport  equations  for  mass,  momentum,  energy,  and  species  transport  are 
presented  here  in  a  Reynolds-Averaged  Navier-Stokes  (RANS)  framework. 

3.1  Navier-Stokes  Governing  Equations 

Conservation  of  Mass: 


dp_  dpuj_  = 
dt  dxj 


(24) 


Conservation  of  Momentum: 


dpu[  dputUj  _  dP  d  T  fduj  duj\ 

dt  dxj  dxi  dxi  V  Vdxj  dxi) 


(25) 
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The  last  term  is  denoted  as  the  Reynolds  stresses  of  the  system,  =  —pu[Uj.  It 
needs  to  be  modeled  to  provide  mathematical  closure  and  to  account  for  turbulence 
effects.  Turbulent  viscosity  formulation  in  RANS  typically  requires  the  modeled 
Reynolds  stress  for  the  Standard  k  —  e  and  RNG  to  be  represented  as, 

Tij  =  pu[Uj  =  2 ptSij  -  1 8ij  ( pk  +  ■  (26) 

Hence,  turbulent  viscosity  is  defined  as  pt  =  C^pik2 /e),  where  k  and  e  are 
calculated  classically  by  transporting  2  representative  equations. 

3.2  Spray  CFD  Models 

Spray  modeling  is  treated  using  the  “blob”  injection  method  of  Bravo  and  Kweon 
(2014a).  Blobs  of  a  characteristic  size  are  injected  following  a  statistical 
distribution  into  the  computational  domain.  Primary  and  secondary  breakups  are 
subsequently  simulated  based  on  the  Kelvin-Helmholts  (KH)  and  Rayleigh-Taylor 
(RT)  instability  methods.  Note  that  the  breakup  length  is  not  determined  a  priori 
(breakup-length  concept)  and  is  calculated  as  a  part  of  the  solution. 

In  the  KH  wave  model,  atomization  is  treated  using  stability  analysis  for  liquid  jets. 
The  breakup  of  the  injected  blobs  and  resulting  drops  of  radius  r0  are  calculated  by 
assuming  that  the  drop  radius  is  proportional  to  the  wavelength  of  the 
fastest-growing  unstable  surface  wave  AKH. 

It  is  written  as 


r  —  B0Akh 


(27) 


where  B0  is  a  model  constant.  The  droplet  size,  and  its  change  of  radius,  is  in  the 
following  way, 


drp  _  _  (r0-r) 

dt  T  kh 


(28) 


where  the  breakup-time  constant,  tkh  ,  is  calculated  as 


3.726£1r0 


(29) 


and  the  maximum  growth  rates  Hkh  and  corresponding  wavelengths  AKH  have 
been  simplified  and  defined  as  follows: 


^KH 


0.34+0. 38We^'5 
(l+Z)(l+1.4r0-6) 


(30) 
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and 


Akh  =  (l+0.45Z05)(l+0.4T°-7) 

a  '  (l+0.87We^,67)°’6 


(31) 


where 

Z  =  We?  5/Rel ,  T  =  We°  5,  Wet  =  ptU2a/a,  Weg  =  pgU2a/a,  and  Ret  =  Ua/vt. 
The  present  RT  mechanism  formulation  includes  viscosity  variations  in  the  growth 
rate  equation, 


mrt  —  ~ 


k2 

K, 


RT 


kRT(^)a.]^  +  kt  f^A2 

\Pl+Pg/J  \Pl+Pg)  Pl+Pg  KI  \Pl+Pg) 


(32) 


where  kRT  is  the  wave  number,  pt  is  the  liquid  viscosity,  pg  is  the  gas  viscosity, 
Pi  is  the  liquid  density,  pg  the  gas  density,  a  is  the  deceleration  of  the  drop,  and  a  is 
the  liquid  surface  tension. 


4.  Computational  Results 

In  Section  4.1  the  2  models  are  systematically  demonstrated  and  compared.  The 
ROM  is  also  compared  to  a  transient  CFD  FOM  described  in  (Senecal  et  al.  2007). 
Section  4.2  compares  the  ROM  output  against  laboratory  measurements  using  a 
mass  threshold  of  90%  that  would  correspond  to  the  penetration  depth  at  which 
90%  of  the  fuel  mass  remains  in  the  liquid  phase.  The  steady  state  also  fluctuates 
with  time,  so  to  compare  with  ROM  output  the  steady-state  measurement  and  FOM 
results  are  averaged  time  or  ensemble  averaged. 

4.1  ROM  and  FOM  Comparison 

The  FOM  simulations  discussed  in  this  section  are  all  modeled  following 
experimental  specifications  for  constant-volume  spray  vessels.  To  get  a  baseline 
before  comparing  to  the  FOM  established  in  Section  3,  we  first  compare  the  ROM 
to  output  from  MoSES,  a  FOM  described  in  (Senecal  et  al.  2007).  Figures  3  and  4 
show  how  the  steady-state  and  transient  parts  of  a  ROM  for 
n-dodecane  compare  against  the  referred  study.  Good  agreement  is  obtained 
between  the  simulations  and  measurements,  so  Fig.  3  uses  the  validation  data  as  a 
proxy  for  the  (Senecal  et  al.  2007)  FOM.  For  this  case,  the  ROM  is  run  for  a  nozzle 
diameter  of  0.246  mm,  a  chamber  temperature  of  1000  K,  and  an  injection  pressure 
of  138  MPa. 


Approved  for  public  release;  distribution  is  unlimited. 

12 


Chamber  Density  (kg/m3) 


Fig.  3  Comparison  of  steady-state  ROM  with  FOM  described  by  Senecal  et  al.  (2007); 
FOM  output  is  taken  at  the  validation  data 


Time  (ms) 


Fig.  4  Comparison  of  transient  ROM  with  FOM  and  measurements  described  by  Senecal 
et  al.  (2007)  for  chamber  density  of  14.8  kg/m3;  FOM  output  is  taken  at  98%  mass  threshold, 
and  measured  value  represents  a  steady  state  of  observation 

For  each  steady  state,  the  results  are  within  a  few  percentage  points  of  each  other. 
Likewise,  the  transient  behavior  of  both  models  generally  agrees,  though  the 
(Senecal  et  al.  2007)  FOM  appears  to  converge  faster  than  the  ROM. 
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Moving  now  to  the  FOM  presented  in  this  study,  5  HPC  simulations  were  executed 
on  the  HPC  to  compare  with  ROM  output  at  steady  state.  For  comparison  to  the 
ROM,  chamber  pressure  was  varied  over  the  5  cases,  which  were  otherwise 
identical.  The  geometry  used  was  a  constant  volume  spray  vessel.  Its  dimensions 
and  conditions  were  as  follows: 


Chamber  diameter 
Chamber  length 
Nozzle  diameter 
Chamber  gas 
Fuel 

Initial  ambient  temperature 
Initial  fuel  temperature 


110  cm 
110  cm 
0.15  mm 
Nitrogen  (N2) 
n-Dodecane 
373  K 
850  K 


Figure  5  shows  the  fuel-injection  process  through  cut  planes  of  temperature 
contours  for  various  increasing  times  (left  to  right)  and  decreasing  chamber 
pressures  (top  to  bottom).  The  model  captures  faster  and  longer  penetration  lengths 
with  lower  chamber  pressures  while  also  showing  higher  turbulence  fluctuations  at 
the  same  conditions. 

0.5  ms  1.0  ms  1.5  ms  2.0  ms  2.5  ms 

1450  kPa 

4450  kPa 


6450  kPa 


9450  kPa 


12450  kPa 


Kelvin:  600  640  680  720  860  800  840  880 

Fig.  5  Transient  behavior  of  n-dodecane  spray  simulation  as  shown  through  temperature 
contours;  liquid  fuel  is  represented  in  black 
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Table  5  shows  the  liquid  penetration  sensitivity  to  fuel-mass  threshold  and  results 
of  5  simulations  with  the  previously  described  conditions. 

Table  5  Simulated  liquid  lengths  (mm)  for  different  mass  thresholds 


Pressure 

(kPa) 

90% 

95% 

97% 

99% 

1450 

22.46 

25.82 

28.11 

32.57 

4450 

15.76 

17.90 

19.21 

21.50 

6450 

15.77 

17.81 

19.02 

21.00 

9450 

12.95 

14.42 

15.30 

16.76 

12450 

10.76 

12.04 

12.83 

14.18 

Figure  6  displays  the  results  of  the  steady-state  ROM  against  the  FOM  outputs  at 
corresponding  densities.  The  best  agreement  between  the  2  occurs  at  the  90%  mass 
threshold,  except  at  the  lower  densities.  The  FOM  also  displays  an  unexpected  dip 
at  17.63  kg/m3.  If  this  point  had  matched  the  trend  of  the  other  points,  it  seems 
likely  the  ROM  would  cross  it  near  the  90%  threshold  like  the  higher  density  runs. 
As  it  stands,  the  ROM  performs  within  a  range  of  reasonable  outputs  from  the 
FOM,  but  they  do  not  compare  as  consistently  at  the  lower  densities  as  they  do  at 
high  density. 


Fig.  6  Comparison  of  steady-state  ROM  with  FOM  at  several  mass  thresholds 
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4.1.1  ROM  and  Experiment  Comparison:  Sandia  National  Laboratory 
(SNL)  Experiments 

Further  zero-dimensional  simulations  were  carried  out  to  demonstrate  its 
applicability  over  a  wide  range  of  conditions  that  have  been  experimentally 
reported.  The  conditions  include  an  evaporating  spray  with  sensitivities  in  the 
chamber  operating  conditions  from  700  to  1300  K  and  across  various  densities  at 
3.6-59  kg/m3. 

Figure  7,  shows  the  results  from  the  single-component  cetane  fuel  case  named 
“ARL  Model”  (for  the  US  Army  Research  Laboratory)  and  compared  to  the 
database  from  the  Department  of  Energy  Engine  Combustion  Network  taken  at 
Sandia  National  Laboratory.  Results  showed  good  agreement  between  model  and 
experiment  with  less  than  5%  discrepancies  at  lower  chamber  densities.  This  can 
be  discerned  from  both  Fig.  7.  Further,  Fig.  8  demonstrates  the  fuel-temperature 
dependency  is  not  significant  for  high-density  cases  of  30  kg/m3. 


E 

E 


-D 

'3 

cr 


Cetane  Fuel  0  p  3.6  kg/m3  ■ 

Lines  ARL  Model  D  073ka/m3 

Symbols:  SNL  Experiment  p '  y  , 

O  p  14.8  kg/m3, 

V  V  p  30.2  kg/m3 

lk  A  p 59.0 kg/m3' 

iquid  Length  (mm) 

-Pfc  CD  00  O  N) 

O  O  O  O  O 

■  Cetane  Fuel  ^  Ta-700K 

\  Lines  :  ARL  Model  °  Ta=850K 

\  Symbols:  SNL  Experiment  O  Ta=1000K 

\  V  Ta=1150K 

\  A  Ta=1300K 

—1 

20 

0 

Chamber  Temperature  (K) 


Chamber  Gas  Density  (kg/m3) 


Fig.  7  Liquid-length  variations  via  (left)  chamber  temperature,  (right)  chamber  gas  density 
space 
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Fuel  Temperature  (K) 

Fig.  8  Dependency  of  liquid  length  with  fuel  temperature 


4.1.2  Impact  of  Fuel  Properties 

Similar  analysis  was  performed  with  the  various  fuel  library  databases  developed. 
In  this  case,  we  utilized  the  compositions  for  cetane,  dodecane,  and  tetradecane  to 
demonstrate  the  fuel  effects  on  the  penetration  parameters.  Figure  9  shows  cetane 
provides  the  best  prediction  of  the  spray  behavior  consistently  over  the  operating 
envelope.  However,  all  fuels  are  able  to  capture  the  behavior  and  dependencies  on 
chamber  temperature.  Further,  we  investigated  the  dependencies  on  the  normalized 
liquid  length  and  B  number  parameters  and  we  were  able  to  demonstrate  a  similarity 
law  that  collapses  all  fuels  under  a  single  curve  as  shown  in  Fig.  10. 
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Symbols:  SNL  experiment  DF2 
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Fig.  9  Fuel  effects  in  liquid  penetration  lengths  over  a  wide  operating  envelope 
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Fig.  10  Similarity  in  normalized  liquid  length  and  B  number  physical  properties 


4.2  ROM  versus  Experimental  Comparison  (ARL  Experiments) 

Figures  1 1  and  12  show  the  scaling  behavior  of  liquid  length  with  3  injector  nozzles 
for  40-,  100-,  and  147-pm  diameter  range.  Figure  11  shows  the  fuel  effects  for 
cetane,  dodecane,  and  tetradecane  versus  chamber  density  for  the  3  injectors.  Figure 
12  shows  the  results  taking  an  optimized  JP-8  surrogate  consisting  of  82% 
dodecane  and  18%  tetradecane. 


Single  Component  Fuels 


Symbols:  Measurement  * 

Lines:  Siebers  Model  v  a-^UUrn 


°0  10  20  30  40  50  60 

Chamber  Density  [kg/m3] 


Fig.  11  Nozzle  diameter  and  fuel  effects  in  liquid  penetration  lengths  over  a  wide  operating 
envelope 
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JP8:  82%dodecane,1 8%tetradecane 


Symbols:  Measurement 
Lines:  Siebers  Model 


0  d=40um 
□  d-IOOum 
O  d=147um 

black:  surrogate 


Chamber  Density  [kg/rrr] 


Fig.  12  Nozzle  diameter  with  JP-8  surrogate  in  liquid  penetration  lengths  over  a  wide 
operating  envelope 

Figure  13  shows  the  log-normal  scale  behavior  between  experiments  and  the  model 
with  good  agreement  for  the  JP-8  surrogate  mixture. 


JP8:  82%dodecane,1 8%tetradecane 


Fig.  13  Log-normal  behavior  of  JP-8  spray 
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4.3  ROM  versus  Real  Engine  Measurements 


The  work  of  Espey  and  Dec  (1995)  provides  a  case  to  validate  the  ROM  against  a 
real  diesel  engine  in  the  form  of  a  direct-injection  Cummins  N-14  with  a  Cummins 
CELECT  injector.  Table  6  shows  the  basic  inputs  needed  to  compare  to  the  ROM. 


Table  6  Basic  inputs  for  Cummins  N-14  diesel  engine 


Nozzle  Diameter 

0.15  mm 

Engine  speed 

1200  rpm 

Injection  pressure 

68  MPa 

Fuel 

67.6%  hepty lmethy lnonane 

32.4%  n-cetane 

To  compare  with  experimental  data,  the  ROM  was  run  for  the  TDC  chamber 
conditions.  An  equal  blend  of  n-dodecane,  n-tetradecane,  and  n-cetane  were  used 
as  a  surrogate  fuel  for  the  mixture  studied  in  experiment.  Figure  14  shows  the 
steady-state  liquid  length  over  12  cycles  for  14  engine  configurations,*  while 
Fig.  15  shows  how  average  transient  behavior  compares. 


Fig,  14  ROM  vs.  real  engine  observed  at  steady-state:  density  on  left  is  16.6  kg/m3, 
temperature  on  right  is  992  K;  error  bars  show  observed  range  of  values  over  12  cycles,  while 
central  points  mark  average  of  all  cycles  and  crank  angles  for  given  TDC  condition 


*  One  of  the  configurations  (pa  =  16.6  kg/m3,  Ta  =  992  K)  is  present  in  both  data  sets. 
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Crank  Angle  (ATDC  )  Crank  Angle  (ATDC  ) 

Fig.  15  ROM  vs.  real  engine  during  injection:  density  on  left  is  16.6  kg/m3,  temperature  on 
right  is  992  K;  solid  curves  are  ROM  output  and  dotted  curves  represent  12-cycle  average  of 
(Espey  and  Dec  1995)  data  at  each  crank  angle 

Each  steady-state  result  from  the  ROM  is  within  the  range  of  observed  values. 
Some  discrepancy  should  be  expected,  as  the  ROM  is  based  on  experimental 
conditions  measured  at  TDC,  while  liquid  lengths  are  measured  several  degrees 
before  TDC.  This  means  the  ROM  is  based  on  chamber  conditions  that  slightly 
overstate  the  density  at  injection.  Were  the  true  injection  conditions  known  the 
ROM  results  would,  therefore,  be  somewhat  lower  than  those  calculated.  Even  so, 
the  ROM  results  are  within  a  few  percentage  points  of  the  observed  average  in  each 
case.  The  temperature  dependence  is  notable  because  while  the  constant-density 
measurements  show  decreasing  length  with  increasing  temperature,  they  do  not 
exhibit  the  smooth,  upward  second-derivative  characteristic  of  the  ROM.  The 
density  dependence  does  exhibit  this  behavior  and  largely  mirrors  that  of  the  ROM, 
though  with  a  slightly  lower  slope. 

The  results  of  Espey  and  Dec  (1995)  are  also  given  by  crank  angle,  allowing  for  a 
comparison  of  transient  behavior,  as  given  in  Fig.  15.  Comparing  the  2  charts 
demonstrates  a  feature  of  the  penetration  correlation  implemented  in  the  ROM. 
Because  Eq.  23  is  dependent  only  on  density  and  nozzle  size,  the  ROM  exhibits  no 
transient  dependence  on  temperature.  This  is  consistent  with  the  measurements, 
which  also  show  temperature  independence,  at  least  until  steady  state.  The  ROM 
follows  transient  measurements  reasonably,  though  experiment  shows  a  more  linear 
trend  for  the  first  moments  after  injection.  The  ROM  loses  its  linearity  much  sooner, 
suggesting  some  limits  to  the  nonvaporizing  spray  model. 
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5.  Discussion 


Compared  to  FOMs,  the  ROM  requires  relatively  few  inputs  to  run,  divided  into  3 
categories.  First  are  settings  for  the  engine,  such  as  engine  speed,  nozzle  design, 
and  injection  parameters.  The  second  category  defines  the  composition  of  fuel  and 
air  components.  The  last  category  defines  chamber-operating  conditions,  which 
include  a  time  variable  and  the  state  variables  of  the  fuel  and  air.  Several  other 
parameters  are  adjustable,  such  as  the  error  thresholds.  These,  however,  are  more 
global  constants  that  do  not  depend  on  the  engine  conditions  being  studied. 

Setup  for  the  FOM  is  comparatively  quite  complex.  Beyond  the  need  for  HPC,  full¬ 
dimensional  CFD  requires  numerous  input  files  to  specify  parameters  such  as 
engine  geometry,  spray  specifications,  and  model  numerical  parameters.  This 
allows  for  high  customization  but  requires  much  more  input  from  the  user.  If  all 
that  is  needed  is  an  estimate  of  liquid  penetration,  executing  a  CFD  model  may  be 
excessive  when  a  simplified  model  will  reproduce  many  of  the  same  results. 

Reasonable  agreement  was  found  between  the  higher-  and  lower-order  models  for 
the  conditions  studied  in  this  report.  The  results  of  the  Espey  and  Dec  (1995) 
measurements  are  also  encouraging  since  they  display  good  agreement  even  during 
combustion,  at  least  until  combustion  products  begin  to  obscure  the  spray.  These 
results  are  of  particular  value  to  test  the  transient  model,  since  there  may  be  concern 
that  a  nonvaporizing  spray  model  will  not  adequately  capture  the  observed 
behavior.  Indeed,  to  approximate  liquid  behavior,  the  velocity  coefficient  Cv  must 
be  modified  to  around  0.22.  Preliminary  agreement  with  measurement  and 
full-order  simulation  suggests  this  method  is  a  reasonable  first  approximation  until 
the  penetration  correlation  can  be  adjusted  for  vaporization. 

To  examine  how  the  initial  conditions  affect  the  ROM,  note  that  it  estimates  the 
steady- state  liquid  length  from  the  product  of  2  factors.  The  first  is  a  normalized 
term 


(33) 


The  second  is  the  length  scale  defined  in  Eq.  19: 


(34) 
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Other  than  a  simple  linear  dependence  on  nozzle  size,  the  scaling  law  is  most 
affected  by  the  temperature,  pressure,  and  density  of  the  chamber.*  The  next  few 
figures  demonstrate  the  effect  of  different  chamber  conditions  on  the  ROM  model 
outputs.  Figure  16  shows  the  effects  gas  density  and  temperature  on  calculations  of 
the  B  number  (evaporation  constant).  Figure  17  shows  how  those  values  translate 
into  a  normalized  length. 


Chamber  Density  (kg/rrr ) 


Fig.  16  Evaporation  constants  for  cetane  injected  into  a  nitrogen  chamber 


Fig.  17  Normalized  liquid  lengths  for  cetane  injected  into  a  nitrogen  chamber 

The  fuels  used  will  also  affect  the  ROM’s  behavior.  Figure  18  shows  how  the 
evaporation  constant  reacts  to  different  pure  fuels  and  also  an  equal  blending  of 
them.  In  this  case  the  blend  is  almost  indistinguishable  from  tetradecane  in  the 
evaporation  constant.  But,  as  Fig.  19  shows,  when  penetration  is  compared,  the 
blend  separates  from  tetradecane  into  a  curve  that  is  still  between  the  fuels. 


Density  is  implicit  in  the  factor  of  tan  a/2,  which  is  computed  by  an  empirical  function  in  Eq.  17. 
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Chamber  Temperature  (K) 


Chamber  Density  (kg/m3) 


Fig.  18  Evaporation  constants  for  3  pure  fuels  and  1  equal  blend;  in  both  cases  the  blended 
evaporation  constant  is  a  simple  average  of  that  for  the  3  pure  fuels — and  since  tetradecane  is 
close  to  average  in  each  case,  a  significant  overlap  results 


Chamber  Temperature  (K)  Chamber  Density  (kg/m3) 


Fig.  19  Steady-state  liquid  lengths  for  3  pure  fuels  and  1  equal  blend 

While  the  ROM  performs  comparably  to  the  cases  studied  in  this  report,  the  ROM 
does  have  limitations.  An  engine  model  studied  by  Scarcelli  et  al.  (2016)  provides 
benchmark  experimental  measures  for  a  gasoline  direct-injection  case  study. 
Scarcelli  et  al.  found  consistent  results  from  the  same  FOM  described  in  Section  3. 
The  FOM  behaves  reasonably,  with  liquid  penetrating  smoothly  up  to  130  mm  into 
the  cylinder  on  each  cycle.  However,  it  was  observed  through  various 
implementations  with  our  models,  that  the  ROM  performs  poorly  under  the 
gasoline  engine  operating  conditions  (over-predicts  by  a  factor  of  more  than  10). 
The  ROM  would  therefore  not  be  recommended  for  such  gasoline  direct  injection 
engines  as  it  stands  in  this  report.  We  suspect  the  low  densities  and  pressures  found 
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during  gasoline  injection*  are  outside  the  range  of  the  empirical  functions  used  in 
Section  2.2.  Future  work  could  explore  this  discrepancy  further  and  modify  the 
model  should  gasoline  behavior  be  desired. 

Other  improvements  to  the  ROM  might  include  the  effects  of  real  gases.  Though 
some  nonideal  effects  are  already  captured  through  the  use  of  compressibility 
factors,  Kurvers  and  Luijten  have  shown  the  simple  division  of  Eq.  8  does  not 
adequately  account  for  the  mixture  of  fuel  and  air.  The  process  they  have  developed 
could  be  used  to  correct  the  model,  most  prominently  at  higher  densities. 

The  transient  segment  of  the  ROM  could  also  benefit  from  further  grounding  in 
theory.  Even  with  the  modification  to  the  penetration  rate,  the  Naber-Siebers 
correlation  is  only  valid  for  a  nonvaporizing  spray  under  the  conditions  noted  by 
Naber  and  Siebers  (1996).  That  the  liquid  fuel  vaporizes  is  a  complication  that  so 
far  has  been  accounted  for  only  by  modifying  the  spray  velocity.  That  experiment 
shows  higher  linearity  than  the  ROM,  which  suggests  a  model  with  more  realistic 
physics  may  provide  better  agreement  with  the  data. 

6.  Conclusions 


The  objective  of  this  research  was  to  establish  a  novel  ROM  for  engine  liquid- 
length  scaling  analysis  and  to  assess  its  validity  over  ranges  relevant  to  Army 
vehicle-propulsion  platforms.  The  formulation  stems  from  the  well-known  Siebers 
model  for  diesel  sprays  that  applies  heat  and  mass  transfer  principles  for  zero¬ 
dimensional,  mixing-limited  conditions.  This  work  extends  this  foundational  model 
to  include  1-D  transient  framework  with  multiphysics  capability.  Fuel  libraries 
have  been  developed  including  pure  fuels  like  n-dodecane,  cetane,  and  tetradecane 
and  also  surrogate  mixtures  to  emulate  jet-propellant  (i.e.,  JP-8)  fuel  properties. 
Companion  CFD  simulations  resolving  the  transient  3-D  spray  behavior  have  been 
performed  to  validate  the  model.  Further  vetting  was  conducted  against 
measurements  for  the  various  cases  of  interest.  The  cases  include  an  evaporating 
single-plume  spray  and  a  single-cylinder  moving  piston  case  near  top  dead  center 
at  diesel-engine  conditions.  The  ROM  provides  a  real-time  engineering  analytical 
tool  for  liquid-length  scaling  that  may  be  used  toward  optimizing  engine 
performance. 

In  summary,  this  report  has  covered  the  following  research: 


*  Diesel  injection  occurs  during  compression  near  TDC  while  gasoline  injection  occurs  during  expansion 
closer  to  bottom  dead  center.  As  a  result,  the  diesel  sprays  studied  tend  to  have  pressures  on  the  order  of  10 
bar  while  the  gasoline  spray  showed  typical  pressures  on  the  order  of  0.1  bar,  a  factor  of  100  lower. 
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.  A  1  -D,  engine-scaling  ROM  was  implemented  based  on  the  work  of  Siebers 
and  Naber  and  significantly  enhanced  to  include  multifuel  effects  and 
transient  injection  physics. 

.  An  FOM,  based  on  a  3-D  CFD  model,  for  diesel  fuel-injection  spray  was 
developed,  implemented,  and  executed  on  the  Army  Engineer  Research  and 
Development  Center’s  Garnet  HPC  system. 

•  The  ROM  was  validated  with  reference  data  for  2  conditions:  diesel  fuel 
injection  in  a  constant  volume  chamber  and  highly  transient  moving  piston 
case. 

.  Model  sensitivity  studies  were  performed  based  on  varying  the  combustor 
operating  conditions  (temperature,  pressure,  and  density),  nozzle  fuel 
temperature,  and  fuel-injector  geometry  (orifice  diameter)  and  showed 
accurate  results. 

•  Suggestions  have  been  made  for  future  improvement  such  as  modeling  the 
effect  of  real  gases  at  supercritical  conditions,  multicomponent 
vaporization,  and  more  accurate  thermodynamic  tables  for  extended  fuel 
libraries  and  operating  envelopes  (such  as  gasoline  fuel  injection). 

The  ROM  remains  largely  valid  in  the  ranges  of  interest,  particularly  in  conditions 
common  to  military-operated  diesel  engines,  though  it  becomes  limited  at  the 
lowest  densities.  Despite  its  simplifications,  it  provides  rapid  and  accurate  results 
for  liquid-spray  analysis  where  experiment  and  intense  CFD  computations  would 
be  costly. 
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Appendix  A.  Stability  and  Convergence  Speed 
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As  mentioned  previously  in  the  reduced-order  model  (ROM)  methodology,  the 
algorithm  used  to  find  the  evaporation  constant  is  based  on  an  iterative  binary 
search.  Just  as  a  binary  search  on  an  array  requires  that  the  array  be  ordered,  the 
root-finding  algorithm  used  here  requires  that  B  be  monotonic  in  Ts  so  that  the 
guess  can  unambiguously  be  updated  based  on  the  value  it  returns. 

Figure  A-l  gives  an  idea  of  how  quickly  the  ROM  converges  given  several 
tolerance  thresholds  £  computed  in  this  equation: 


Low  values  for  e  allows  for  higher  confidence  in  B,  since  the  difference  Bg  —  Bh 
must  be  smaller  to  achieve  it.  Around  the  10-3  tolerance  level,  any  further  gains  in 
precision  are  much  smaller  than  inconsistencies  with  experimental  data  and  the 
FOMs,  so  higher  precision  than  that  is  unlikely  to  yield  more  reliable  results. 

The  method  converges  in  logarithmic  time,  matching  what  would  be  expected  from 
a  binary  search.  The  consistency  of  convergence  time  should  allow  for  a  fairly 
accurate  estimate  of  the  time  required  to  run  the  code.* 


O  u 
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Fig.  A-l  Run  time  to  determine  evaporation  constant  is  logarithmic  in  the  tolerance  required; 
higher  precision  is  shown  from  left  to  right  in  orders  of  magnitude,  which  requires  lower 
tolerance 


*  For  context,  this  code  was  timed  on  a  Dell  Latitude  laptop  (Intel  processor  at  2.5  GHz)  to  run  about 
0.5  ms  per  iteration.  A  tolerance  threshold  of  8  =  10-3  or  8  =  10-4  was  used  for  most  of  the  results  in  this 
study,  translating  to  around  12  ms  to  determine  B.  Running  the  rest  of  the  model  takes  on  the  order  of  a  few 
milliseconds.  Many  of  the  condition  sweeps  shown  in  the  rest  of  this  paper  were  created  in  less  than  a  second. 
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Figure  A-2  shows  the  behavior  of  the  evaporation  constant  (or  B  number)  as  a 
function  of  the  saturation  temperature  in  the  ROM  code.  The  model  shows  fuel 
sensitivities  at  the  lower  saturation  temperatures  up  to  475  K. 


Fig.  A-2  Evaporation  constant  dependence  to  saturation  temperature  for  various  single 
component  fuels. 


Figures  A- 3  and  A-4  show  the  behavior  of  the  enforced  physical  functions — right- 
hand  side  (RHS)  and  left-hand  side  (LHS) — used  in  the  ROM  code  formulation 
with  saturation  temperature.  Further,  the  model  error  calculated  as  the  difference 
between  the  functions  and  fuel  sensitivities  is  demonstrated. 


Fig.  A-3  Behavior  of  the  RHS  and  LHS  functions  with  saturation  temperatures  of  various 
single-component  fuels. 
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Fig.  A -4  Model  error  (difference  between  RHS  and  LHS)  with  saturation  temperatures  of 
various  single-component  fuels 
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Appendix  B.  MATLAB  Code  for  Reduced-Order  Model  (ROM) 


This  appendix  appears  in  its  original  form,  without  editorial  change. 
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The  MATLAB  code  is  executed  through  various  modules  and  functions,  which 
contain  the  algorithmic  formulation,  solution  procedure,  and  thermodynamic 
properties. 

-k'k-k-k-k-k-k-k-k'k'k'k'k'k'k 

L  L_ROM_Ru  n  T  e  mp  late 

%%  FLOW  OF  INFORMATION  BETWEEN  FILES 

clear,  clc 

%  Load  inputs 
LL_ROM_Input Values 

%  [MODIFICATIONS  TO  AMBIENT  COMPONENTS  &  TEMPERATURE] 

%  Calculate  ambient  properties  common  to  all  fuels 
AmbientPropert ies 
Ambient Mixture 

%  [MODIFICATIONS  TO  OTHER  CHAMBER  CONDITIONS] 

%  Calculate  liquid  length  and  other  properties  for  each  case 
LL_ROM_B 1 ends olut ion 

%  [MODIFICATIONS  TO  TIME  &  ENGINE  DESIGN] 

%  Calculate  transient  liquid  length 
LL_ROM_Transient Solution 


-k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

LL_ROM_Input Values 

%%  UNIVERSAL  CONSTANTS 

R  =  8.314E-02;  %  Ideal  Gas  Constant  (bar-m3/kmol-K) 

a  =  0.66;  %  Correction  factor  for  real  spray 

b  =  0.41;  %  Correction  factor  for  simplified  scaling  law  assumptions 
%%  ENGINE  SETTINGS 

rpm  =  1200;  %  Engine  speed  (rev/min) 

%  Nozzle 

d  =  150E-06;  %  Nozzle  diameter  (m) 

NozzleAngle  =  0;  %  Nozzle  angle  (degrees  from  cylinder  axis  -  don’t 
confuse  with  spray  angle) 

Ca  =  0.80;  %  Area-contraction  coefficient  (default  is  average  of 

Siebers  1999  values) 

Cv  =  0.22;  %  Velocity  coefficient  (default  is  fit  to  Espey  &  Dec 

1995  data) 
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c  =  0.2640;  %  Spray  model  parameter  (default  is  average  of  Siebers 
1999  values) 

%  Injection 

Pinj  =  680;  %  Fuel  injection  pressure  (bar) 

CAinj  =  -11.25;  %  Crank  angle  at  fuel  injection  (deg  ATCD) 

%%  FUEL  &  AIR  COMPONENTS 


%  Mole  fractions  of  fuel  components: 
Yf  =  [ 


1.0,  .. 

.  "6 

Dodecane : 

fuel_f lag 

0.0, 

.  "6 

Tetradecane : 

f uel_f lag 

0.0  .  . 

] ; 

o. 

Cetane : 

f uel_f lag 

nFuels  =  max (size (Yf ) )  ; 


%  Mole  fractions 
Ya  =  [ 

0 .7192,  .  .  . 

0.0000,  ... 
0.1923,  ... 

0.0885  ... 


of  ambient  air  components: 

%  Nitrogen 
%  Oxygen 

%  Carbon  Dioxide 
%  Water  Vapor 


] ; 


%%  CHAMBER  OPERATING  CONDITIONS 
t  =  0.3E-03;  %  Time  from  injection  (s) 

%  Fuel 

Tf  =  327;  %  Initial  fuel  temperature  (K) 

%  Ambient  air  (make  sure  all  3  values  follow  gas  law) 
Ta  =  700;  %  Initial  ambient  temperature  (K) 

RhoA  =  0.344;  %  Initial  ambient  density  (kg/m3) 

Pa  =  0.458;  %  Initial  ambient  pressure  (bar) 

%%  CALCULATION  OPTIONS 

epsilon  =  0.001;  %  Relative  error  threshold  for  B 


^^^^^itizitizitizitizitititizit'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

'k'k'k'k'k'k-k-k-k-k-k'k-k'k'k 

Ambient Properties 

%  Sets  properties  for  ambient  air  components 
%  From  API  data  (Table  1C3) 

AirProps  =  ... 


MW  , 

Tc  , 

Pc  , 

Vc  , 

Zc  ; 

28 . 02, 

126.2, 

O 

O 

00 

0 . 0891, 

0.289; 

o. 

Nitrogen 

32 . 00, 

154 . 6, 

50.43, 

0 . 0733, 

0.288; 

o. 

Oxygen 

^ — i 

o 

304.2, 

00 

00 

00 

[> 

0.0939, 

0.274; 

...  "6 

Carbon  Dioxide 

CM 

O 

00 
t — 1 

647 .1, 

220.6, 

0 . 0559, 

0.229; 

...  "6 

Water  Vapor 
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[INSERT  EXTRA  COMPONENT  VALUES  AS  NEEDED] 


] ; 


%  Molecular  weights 

MWas  =  AirProps ( : , 1 ) ;  %  kg/kmol 


%  Values 
Tcras  = 
Pcras 
Veras 
Z  c  r  a  s 


at  Critical 
AirProps ( : , 
AirProps ( : , 
AirProps ( : , 
AirProps ( : , 


Point 


2) ; 

"6 

Temperature  (K) 

3)  ; 

"6 

Pressure  (bar) 

4) ; 

"6 

Volume  (m3/kmol) 

5)  ; 

o. 

Compressibility 

clear  AirProps 


'k'k'k'k'k'k-k-k'k'k-k-k-k-k'k 

AmbientMixture 


%%  COMBINED  VALUES  FOR  AMBIENT  MIXTURE 

%  Molecular  Weight 

MWa  =  Ya*MWas ;  %  kg/kmol 


%  Pseudocrit ical 
Zcra  =  Ya*Zcras; 
Tcra  =  Ya^Tcras; 
Vera  =  Ya*Vcras; 
Pcra  =  Ya*Pcras; 
%  Pcra  =  Zcra  * 


values  (Modified  Prausnit z-Gunn  Rule) 
%  Compressibility 
%  Temperature  (K) 

%  Volume  (m3/kmol) 

%  Pressure  (bar) 

R  *  Tcra  /  Vera;  %  bar 


Za  =  1  ; 


'k'k'k'k'k'k-k-k-k-k-k-k-k-k'k 

LL_ROM_B 1 ends olut ion 

%%  LL_ROM_B lends olut ion  CALCULATES  THE  LIQUID  LENGTH  OF  A  BLENDED 
FUEL  SPRAY 

%  Ensure  Yf  &  Ya  normalized: 


if  sum ( Yf ) = 
Yf / sum ( Yf ) ; 

=0; 

end 

error ( 1  No 

fuel  components 

selected ' ) ; 

else  Yf 

if  sum ( Ya) = 
Ya/sum ( Ya) ; 

=0; 

end 

error ( ’ No 

air  components 

selected ’ ) ; 

else  Ya  = 

B  =  0;  RhoF  =  0; 

%  Initial  Enthalpy  of  Chamber 

Ha  =  AmbientEnthalpy (Ta) ;  %  kJ/kmol 

Par  =  Pa/Pcra;  %  Reduced  ambient  pressure 
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global  fuel_flag 

for  fuel_f lag=l : nFuels 

if  Yf ( fuel_f lag)  ==  0;  continue;  end  %  Skips  fuel  if  component 
is  zero 

FuelPropert ies 

%  Initial  Enthalpy  of  Liquid  Fuel 
Hf  =  FuelEnthalpy (Par ) ; 

%%  RUN  PURE  FUEL  SOLUTION  FOR  EACH  FUEL 


%  Search  space  for  Ts 
TsMin  =  Tf;  TsMax  =  Ter; 

LL_ROM_PureSolut ion 

%  Vapor  Coefficient  of  Fuel  Blend.  Mean  Evaporation  Constant 
method 

B  =  B  +  Yf (fuel_flag) *FuelB; 

%  Liquid  density  from  Rackett  equation.  From  API  data  (Procedure 
6A2 . 13) 

Tf r  =  Tf/Tcr; 

invRhoF  =  (R*Tcr )  /  (Pcr*MWf)  *  ZraA(l  +  ( 1-Tf r ) A (2/7 ) ) ; 

RhoF  =  RhoF  +  Yf ( fuel_f lag) /invRhoF;  %  kg/m3 

end 

%%  CALCULATE  THETA  BASED  ON  SIEBERS  CORRELATION 
RhoRatio  =  RhoA/RhoF; 

TanHalf Theta  =  c  *  (RhoRat ioA 0 . 1 9  -  0 . 0043/sqrt (RhoRatio) )  ; 

HalfTheta  =  atand (TanHalf Theta)  ; 

%%  CALCULATE  LIQUID  LENGTH  BASED  ON  SIEBERS  SCALING  LAW 

LLnorm  =  b  *  sqrt((l  +  2/B) A2  -  1); 

LLscale  =  1/a  *  d*sqrt (Ca)  /  ... 

(sqrt (RhoRatio)  *  TanHalf Theta) ; 

%  Liquid  Length  as  measured  vertically  down  axis  of  chamber 
LL  =  LLscale  *  LLnorm  *  cosd (NozzleAngle) ;  %  m 

LL_mm  =  LL*1000;  %  mm 

%%  CALCULATE  TIME  SCALES  FOR  TRANSIENT  SOLUTIONS 
%  Engine  speed 

AngVel  =  rpm  *  360  /  60;  %  Degrees/s 
%  Fuel  velocity 

v2  =  2* (Pinj  -  Pa) /RhoF  *  100000;  %  m2/s2 
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FuelVel  =  Cv  *  sqrt(v2);  %  m/s 
tscale  =  LLscale/FuelVel ; 

%  Time  when  LL(t)  =  LL 

LLtnorm  =  Penet rat ionTime (LLnorm) ;  %  Normalized  (from  injection) 

%  Crank  angle  when  LL(t)  =  LL 

CAeqb  =  CAinj  +  AngVel*t scale*LLtnorm;  %  Degrees 


-k-k'k-k'k'k'k'k'k'k'k'k'k'k'k 

LL_ROM_Pure Solution 

%%  LL_ROM_PureSolut ion  CALCULATES  Ts  OF  A  PURE  FUEL  ITERATIVELY 
WITH  A  BINARY  SEARCH  ALGORITHM 

Ts  =  (TsMax  +  TsMin)  /  2;  %  Initial  guess  for  Ts 

TsRange  =  (TsMax  -  TsMin)  /  2;  %  Initial  size  of  Ts  search  space 

iter  =  0; 

accepted  =  false;  failed  =  false; 

while  -accepted 

%%  START  OF  LOOP 
iter  =  iter  +  1; 
if  iter>30 

failed  =  true; 
continue 

end 

%  Vapor  pressure  of  fuel  at  aturation 
Ps  =  FuelSaturat ionPressure (Ts ) ;  %  bar 

Prs  =  Ps/Pcr; 

Trs  =  Ts/Tcr; 

%  Fuel  enthalpy  at  saturation 

Hfs  =  SaturatedFuelEnthalpy (Prs , Trs ) ;  %  kJ/kg 

%  Ambient  enthalpy  at  saturation 

Has  =  100*Ya*R*AmbientEnthalpies (Ts) /MWa; 

%  Fuel  Compressiblity 

Zf  =  FuelCompressibility (Trs )  ; 

%  Evaluate  Vapor  Coefficient  for  Iterative  Solution  (Ts) 

Bg  =  (Za  *  Ps  *  MWf )  /  (Zf  *  (Pa-Ps)  *  MWa) ;  %  B  value  from 

gas  law  (Increases  with  Ts) 

Bh  =  (Ha  -  Has)  /  (Hfs  -  Hf)  ;  %  B  value  from 

enthalpies  (Decreases  with  Ts) 

FuelB  =  (Bg  +  Bh)  /  2;  %  Mean  of  B 

values  adopted 
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%  Calculate  error 

Error  =  (Bg  -  Bh)  /  FuelB; 

%  Adjust  Ts  for  Error  value 
if  abs (Error)  <  epsilon 
accepted  =  true; 

else 

TsRange  =  TsRange/2; 
if  Error  >  0 

Ts  =  Ts  -  TsRange; 

else 

Ts  =  Ts  +  TsRange; 

end 

end 

end 


'k'k'k'k'k'k-k-k-k-k-k-k-k'k-k 

LL_ROM_T r an sient Solution 

%%  LL_ROM_Transient Solution  CALCULATES  THE  LIQUID  LENGTH  AT  A  GIVEN 
TIME 

%  Time  from  injection 

tnorm  =  t/tscale;  %  Normalized 

if  tnorm<=0 

xnorm  =  0; 

elseif  tnorm  <  LLtnorm 
%  Initial  guesses: 
xdomain  =  LLnorm; 
xnorm  =  xdomain  /  2; 

accepted  =  false; 
while  ~accepted 

tguess  =  Penet rat ionTime (xnorm) ; 

Error  =  1  -  tguess/tnorm; 

if  abs (Error ) <0 . 001 
accepted  =  true; 

else 

xdomain  =  xdomain/2; 
if  Error  <  0 

xnorm  =  xnorm  -  xdomain/2; 

else 

xnorm  =  xnorm  +  xdomain/2; 

end 

end 


end 
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else 

xnorm  =  LLnorm; 

end 

%  Transient  liquid  length 
LLt  =  LLscale  *  xnorm;  %  m 


-k-k'k-k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

-k'k-k'k'k'k'k'k'k'k'k'k'k'k'k 

FuelProperties 

%%  SET  PURE  FUEL  PROPERTIES 
%  From  API  data  (Table  1C1  &  Table  6A2.14) 


global  fuel 


switch  fuel_flag 
case  1 


fuel 

=  'Dodecane 

?  . 
t 

MWf 

=  170.34; 

"6 

Molecular  Weight 

(kg/kmol ) 

Ter 

=  658; 

"6 

Critical 

Temperature  (K) 

Per 

II 

M 

00 

ISO 

"6 

Critical 

Pressure 

(bar) 

Zra 

=  0.2471; 

o. 

Rackett 

Parameter 

case  2 


fuel 

MWf 

Ter 

Per 

Zra 

case  3 


’ Tetradecane ' ; 

198.39;  %  Molecular  Weight  (kg/kmol) 

693;  %  Critical  Temperature  (K) 

15.7;  %  Critical  Pressure  (bar) 

0.2270;  %  Rackett  Parameter 


end 


fuel 

=  'Cetane' 

r 

MWf 

=  226.45; 

"6 

Molecular  Weight 

(kg/kmol ) 

Ter 

=  723; 

o. 

Critical 

Temperature  (K) 

Per 

=  14.0; 

o. 

Critical 

Pressure 

(bar) 

Zra 

=  0.2386; 

"6 

Rackett 

Parameter 

[INSERT  EXTRA  COMPONENT  VALUES  AS  NEEDED] 
otherwise 

error ('Not  a  valid  Fuel') 


'k'k'k'k'k-k-k-k-k-k-k-k-k-k-k'k-k'k'k'k'k-k-k'k-k'k-k-k-k'k-k'k'k'k'k'k'k-k-k'k-k-k-k'k-k'k-k-k'k'k'k'k'k'k 

-k'k-k'k'k'k'k'k'k'k'k'k'k'k'k 

SetHaVals 

%  From  NASA  data  (McBride  et  al  1993,  Table  2) 

%  Valid  between  200K  -  6000K  unless  otherwise  specified 

global  HaValsHiT;  global  HaValsLoT ;  global  HaValsZero; 

HaValsHiT  =  [  . . . 
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R2C1 


R2C2 


R2C3 


R2C4 


R2C5  R3C1 

2 . 95257 62 6E  +  00  1 . 39690057E-03 

11  -4 . 60755321E-15  -9 . 23948645E+02 ; 

3. 6  60  9  60  83E  +  0  0  6 . 5636852  3E-0  4 

11  -1 .29913248E-15  -1 . 21597725E  +  03; 

4 . 63  65  94 93E  +  0  0  2 . 74131 991E-03 

10  -9. 16103468E-15  -4 . 9024 9341E+04 ; 

2 . 67703787E  +  00  2 . 97  31832 9E-03 

11  -4 .26900959E-15  -2 . 98858 938E+04 ; 


-4 . 92  631 6  91E-07 
. . .  Nitrogen 
-1 . 41149485E-07 
. . .  Oxygen 
-9. 9582  8531E-07 
.  .  .  Carbon  Dioxide 
-7 . 7  3769 6 90E-0 7 
. . .  Water  Vapor 


7 .  86010367E- 
2 . 057  97 658E- 
1 . 60373011E- 
9 . 44336689E- 


] ; 


HaValsLoT  =  [  . . . 

%  R3C3  R3C4 

R4C2  R4C3 

3 . 5310 052 8E+00  -1 . 23660987E-04 
09  -1 . 40881235E-12  -1 . 04 697 628E+03; 

3 . 78245636E+00  -2 . 99673415E-03 
09  3 . 24372836E-12  -1 .  0  63 94  35 6E  +  03; 

2 . 35677352E+00  8 . 9845 967 7E-03 

09  -1 . 43699548E-13  -4 . 8371 9697E+04 ; 

4 . 198 6405 6E+00  -2 . 03643410E-03 
09  1 . 771 9781 7E-12  -3 . 02 9372 67E+04 ; 

]  ; 


R3C5  R4C1 

-5 . 02  999437E-07  2 . 4  353  0 6 12E- 

. . .  Nitrogen 

9. 84730200E-06  -9 . 6812 9508E- 
. . .  Oxygen 

-7 . 12  35  62  6  9E-0  6  2 . 4  5  919022E- 

.  .  .  Carbon  Dioxide 

6 . 5204021 1E-0  6  -5 . 48797062E- 
. . .  Water  Vapor 


HaValsZero  =  [  ... 

%  R4C5 

0 . 00000000E  +  00; 
0 . 00000000E  +  00; 
-4 . 73281047E+04; 
-2 . 908481 68E+04 ; 
]  ; 


. . .  Nitrogen 
. . .  Oxygen 
. . .  Carbon  Dioxide 
. . .  Water  Vapor 


'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

Ambient Enthalpies 

%  Enthalpy  of  ambient  air,  up  to  factor  of  R 
%  NASA  data  &  formula  (McBride  et  al  1993) 

function  Har  =  AmbientEnthalpies (T) 

global  HaValsHiT;  global  HaValsLoT;  global  HaValsZero; 

%  Set  coefficients 
if  T>=1 0 0 0 

HaVals  =  HaValsHiT;  %  High-temperature  coefficients 

else 

HaVals  =  HaValsLoT;  %  Low-temperature  coefficients 

end 

%  Set  indeterminates 

vars  =  [T;  TA2/2;  TA3/3;  TA4/4;  TA5/5;  1]; 
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%  Array  of  temperature  components  for  each  species 
Har  =  HaVals*vars  -  HaValsZero;  %  kJ/kmol-R 

end 


^^^^^ititizititizitizitititit'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

FuelEnthalpy 

%  Initial  enthalpy  of  liquid  fuel 
%  Curve  fit  by  Schihl  et  al  from  API  data 

function  Hf  =  FuelEnthalpy (Par ) 

global  fuel 
switch  fuel 

case  'Dodecane' 


Cl 

=  1.277; 

C2 

=  316.3; 

case  ' 

Tetradecane 

Cl 

=  1.1718; 

C2 

=  319.2; 

case  ' 

Cetane ' 

Cl 

=  1.0429; 

C2 

=  320.95; 

end 

%  calculate .  .  . 

Hf  =  Cl *Par  +  C2 ; 

end 


^^^^^itititititititit'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

'k'k'k'k'k'k'k-k'k'k-k-k-k'k'k 

Fuel SaturationP res sure 

%  Vapor  pressure  at  saturation  temperature 
%  From  API  data  (Procedure  5A1.1,  Table  5A1.2) 

function  Ps  =  FuelSaturat ionPressure (Ts ) 

global  fuel 
switch  fuel 

case  'Dodecane' 

A  =  170.24; 

B  =  25990; 

C  =  20.822; 

D  =  2 . 9759E-06; 

E  =  645190; 
case  ' Tetradecane ' 
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A  =  130.78; 

B  =  20072; 

C  =  15.743; 

D  =  2.3531E-06; 
E  =  -1008900; 
case  'Cetane' 

A  =  174.20; 

B  =  28534; 

C  =  21.090; 

D  =  2.5228E-06; 
E  =  88111; 

end 


TsR 


1.8  *  Ts;  %  Temperature  (Rankine) 


Ps  =  exp ( 

A 

-  B/TsR 

-  C*log (TsR) 
+  D*TsRA2 

+  E/TsRA2 
) ;  %  psi 


Ps  =  0.06895  *  Ps;  %  bar 


end 


^^^^^ititizizitititizitititit'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

SaturatedFuelEnthalpy 

%  Curve  fit  by  Schihl  et  al  from  API  data 
function  Hfs  =  SaturatedFuelEnthalpy (Prs,  Trs) 


global  fuel 
switch  fuel 


case  1 

'  Dodecane ' 

A 

= 

1562.5; 

B 

= 

444 . 84; 

C 

= 

1921 .9; 

D 

= 

4996.7; 

E 

= 

1978; 

case  1 

'  Tetradecane 

A 

= 

1683.8; 

B 

= 

469.88; 

C 

= 

3712.2; 

D 

= 

8085.6; 

E 

= 

3229.2; 

case  1 

1  Cetane ' 

A 

= 

1869.7; 

B 

= 

550 . 15; 

C 

= 

929.51; 

Approved  for  public  release;  distribution  is  unlimited. 

43 


end 


D 

E 


=  3520.3; 
=  1283.4; 


if  Prs<0 . 2 

Hfs  =  A*Trs  -  B; 

else 

Hfs  =  -C*Trs A2  +  D*Trs 

end 


%  kJ/kmol 
E;  %  kJ/kmol 


end 


'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

FuelCompres sibil it y 

%  Curve  fit  by  Schihl  et  al  from  API  data 

function  Zf  =  FuelCompressibility (Trs ) 

global  fuel 
switch  fuel 

case  'Dodecane' 

A  =  16.85; 

B  =  36.104; 

C  =  26.425; 

D  =  7.5406; 

case  ' Tetradecane ' 

A  =  17 . 924; 

B  =  36.143; 

C  =  24.71; 

D  =  6.6857; 

case  ’Cetane’ 

A  =  16.587; 

B  =  34.594; 

C  =  24.531; 

D  =  6.869; 

end 

%  calculate .  .  . 

Zf  =  -A*TrsA3  +  B*TrsA2  -  C*Trs  +  D; 
end 


'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

PenetrationTime 

%  Normalized  time  required  for  liquid  length  to  reach  normalized 
distance  argument 
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From  Naber  &  Siebers 


function  tnorm  =  Penet rat ionTime (xnorm) 
tnorm 

xnorm/ 2 

+  xnorm/4  *  sqrt(l  +  16*xnormA2) 

+  log(4*xnorm  +  sqrt(l  +  1 6*xnormA2 ) ) / 1 6 ; 


end 
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List  of  Symbols,  Abbreviations,  and  Acronyms 


3-D 

3-dimensional 

API 

American  Petroleum  Institute 

ARL 

US  Army  Research  Laboratory 

CFD 

computational  fluid  dynamics 

FOM 

full-order  model 

HPC 

high-performance  computing 

LHS 

left-hand  side 

KH 

Kelvin-Helmholts 

MEC 

mean  evaporation  constant 

NASA 

National  Aeronautics  and  Space  Administration 

RANS 

Reynolds- Averaged  Navier-Stokes 

RHS 

right-hand  side 

RNG 

renormalization  group 

ROM 

reduced-order  model 

RT 

Rayleigh-T  ay  lor 

SNL 

Sandia  National  Laboratory 

TDC 

top  dead  center 
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